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(57) Abstract: Methods 
and systems are disclosed 
for processing data used for 
hydrocarbon extraction from the 
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a set of predicted critical points. 
The subdivisions are preferably 
generated such that a shape family 
is selected according to a best fit 
between the critical points from 
the sampled data and the predicted 
critical points of the selected shape 
family. 
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System and Method for Representing and Processing and 
Modeling Subterranean Surfaces 

FIELD OF THE INVENTION : 
5 The present invention relates to the field of earth 

models for subterranean surfaces. In particular, the 
invention* relates to systems and methods for improved 
representations and processing techniques for 
subterranean earth surfaces in earth models used in the 
10 exploration and production of hydrocarbon reservoirs. 

BACKGROUND OF THE INVENTION ; 

In the field of processing earth model data for use 
in the extraction of hydrocarbons from the earth, 

15 significant resources have been invested in creating 
functionality . and stabilizing the software quality of 
modeling technology. The efforts have been based on 
faceted representations, and in particular triangulated 
surface methods. However, there are number of key 

20 problems associated with extending triangulated surface 
methods . 

1. Sampling a smooth surface discretely, for 
example with points arranged in a triangle mesh, is 
inherently inefficient. In contrast, smooth surfaces can 

25 be represented as a Taylor series or as an eigen- function 
expansion, e.g., Fourier series of some form. An eigen- 
function expansion can be used to compute an algebraic 
expression to evaluate normal fields, tangent fields, 
etc. This is inherently more compact and efficient than a 

30 2D or 3D sampling lattice with some sort of 
interpolation . 
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2 . The lack of differentiation makes calculation 
of a triangulated surface intersection algorithm 
numerically delicate. 

3.. Triangles can be used to define shape, but 
5 triangles do not efficiently convey an intuitive sense of 
shape. High resolution and high sampling density make the 
problem more difficult. 

4. Low efficiency triangulation forces application 
developers to either worry about memory management or to 

10 curb flexibility of data set processing.* 

5. Efficient management of large models means 
localizing change. It is possible in principle to develop 
localization methods uging triangulated surface models 
but there are 'numerical stability issues seen in 

15 reference to intersection. 

SUMMARY OF THE INVENTION : 

Thus, it is an object of the present invention to 

provide an improved system and method for processing data 
20 • used for .hydrocarbon extraction. 

Advantageously the invention allows for improved 

memory and CPU efficiency of implicit surface 

construction and editing algorithms. 

According to the invention, a method is provided for 
25 processing data used for hydrocarbon extraction from the 

earth. The method includes the following steps. 

Receiving sampled data representing earth structures. 

Identifying one more symmetry transformation groups from 

the sampled data. Identifying . a. set of critical points. 
30 from the sampled data. Generating a plurality of 

subdi visions of shapes, the subdi visions together 

representing the earth structures, the generation being 

based at least in part on the set of identified critical 
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points and the symmetry transformation groups. Processing 
earth model data using the generated subdivision of 
shapes. And, altering activity relating to extraction of 
hydrocarbons from a hydrocarbon reservoir based on the 
5 processed earth model data. 

The identified symmetry transformation group is 
preferably a set of dif f eomorphisms that act on a 
topologically closed and bounded region in space-time 
such that under transformation the region occupies the 
10 . same points in space. 

The identified symmetry transformation groups 
preferably correspond to a plurality of shape families, 
each of which includes a set of predicted critical 
points. The subdivisions are preferably generated such 
15 that a shape family is selected from the plurality of 
shape families that corresponds to the identified 
symmetry transformation group. The selection is 
preferably being based on closeness of correspondence 
* between the identified critical points from the sampled 
20 data and the predicted critical points of the selected 
shape family. 

Each shape family preferably has an associated set 
of symmetry transformation group orbits, with each orbit 
being associated with orbit information that specifies 
25 whether the* orbit contains a predicted critical point and 
value of the Gaussian curvature of a point in the orbit. 
The orbit information from the set of symmetry 
transformation group orbits associated with the selected 
shape farriily is preferably applied to the sampled data 
30 thereby generating a unique specification of a shape from 
the selected shape family. Each of the plurality of 
subdivisions of shapes is preferably generated by 
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identifying a part of the uniquely specified shape that 
corresponds to the sampled data. The identified parts ' 
. are assembled, thereby generating a representation of the 
earth structures . 
5 The plurality of subdivisions are preferably 

. generated such that the number of parameters in each 
subdivision times the number of subdivisions is 
substantially less then would be needed using a faceted 
representation method, and the plurality of subdivisions 
10 is more numerically stable than third order or higher 
representation. 

The invention is also embodied in a system for 
improved extraction of hydrocarbons from the earth, and 
in a computer readable medium capable of causing a 
15 computer system to carry out steps for processing data 
relating to earth structures for the extraction of 
hydrocarbons . 

■ BRIEF DESCRIPTION OF THE DRAWINGS : 

20 

Figure 1 shows an open surface and its embedding in 
a closed surface; 

Figure 2 shows a salt mass's steep flanks and 
overhangs ; 

25 Figure 3 shows an example of a 4D representation of 

a field in Turkmanistan; 

Figure 4 is an image of the MacKenzie River Delta; 
Figure 5 shows some combinations of involved 
spherical harmonic polynomials, presented in spherical 
30 polar coordinates; 

Figure 6 illustrates that in most cases a surface 
evolves under mcf to a point; " 
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Figures 7a- f forms a series of six images showing 
vpmcf suppressing noise; 

Figure 8 illustrates the orthogonality condition of 
the theorem proposed by Athanassenas ; 
.5 Figure 9 is an aerial image of part of Big Bend 

National Park, showing the approximation of a plateau to 
a characteristic length scale cone; 

Figure 10 is a side view of a noise cone structure; 

Figure 11 is a satellite image of the Labrador 
10 Trough; 

Figure 12 shows s a sequence of folded sediment on 
the coast of the Gulf of Oman; 

Figure 13 is an image illustrating progressive 
flattening of an overburden covering a large salt 
15 intrusion; 

Figure 14 is a diagram illustrating the Morse 
theoretical cell decomposition for a simple configuration 
of a capped and bent, cylinder; 

Figures 15-17 show diagrams to aid in the 
20 understanding of the bulls eye construction; 

Figure 18 shows . two views of an example of a monkey 
saddle; 

Figure 19a shows the Reeb graph of a standard torus; 

Figure 19b schematically illustrates a 2D cell 
25 suspension that is* induced from the axes and planes of 
symmetry and critical point theory; 

Figure 20 is cross section of the torus shown in 
Figure 19a; 

Figure 21 is a schematic of the shape synopsis 
30 diagram of the torus shown in Figure 19a; 

Figure 22 is a cyclide that is shaded according to 
Gaussian and mean curvature; 
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Figure 23a is a bi-torus with its associated REEB 
diagram; 

Figure 23b is the visual representation of the bi-. 
torus shape synopsis diagram; 
5 Figure 24 is a diagram of the octree with a- coarse 

level and leaf level shape index relationship indicated; 

The Figure 25 is a diagram shows part of the French 
model ; 

Figure 2 6 is an image of the topography of Crater 
10 Lake, Oregon; 

Figure 27 shows a salt weld in the Gulf of Mexico; 
Figures 28a-c illustrate an example of the misfit 
reduction process ; 

Figures 29a-c illustrates an example of blending a 
15 non-dif f erentiable join of two collars; 

Figure 30 shows as a geological example a water 
breach as indicated by the white arrow; 

Figure 3*1 is a NASA Shuttle Mission photograph of 
the Richat Structure in Mauritania; 
20 Figure 32 shows the natural analog to a conformal 

grid with a proportionally spaced correlation scheme; 

Figure 33 illustrates a non-conf ormal 3D Cartesian 

grid; 

Figure 34 is an image of the Devil's Potholes, South 
25 Africa; 

Figure 35 is an image of the Yukon River delta; 

Figures 3 6a and 36b show, for reference, the 
background Zechstein Salt and the region in the Zechstein 
where the vsp was acquired; - 
30 Figure 37 show the frame graph that ties the vsp 

region of interest to the Zechstein Salt background; 
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Figures 38a, 38b, 39a and 39b illustrate the 
separation of faulted sediments from unfaulted sediments; 

Figure 40 illustrates a time-lapse seismic 
evolution; 

5 . Figure 41 shows the reference structure for spatial 

frames to define a topology graph; 

Figure 42 is a schematic illustration of the 
• definition of the variables in the mean curvature 
estimator; 

10 Figure .43 is an image of the regularly sampled input 

surface representing the present example of the top of 
the salt, rendered by drawing a subset of evenly spaced 
inlines and crosslines; 

Figure 44 shows the, surf ace when re- sampled by 
15 interpolating missing sample points; 

. Figure 45 shows the re-sampled surface 'after 25 
iterations of smoothing; 

Figure 46 shows the shape index map for the re- 
sampled surface after 25 iterations of smoothing; 

2 0 • Figure 47 shows. that there are 33 shapes in the 

example shown in Figure 46; 

Figure 48 is a schematic illustration of a system 
for improved extraction of hydrocarbons from the earth, 
according a preferred embodiment of the invention ; 
25 Figure 49 shows further detail of a data processor 

according to preferred embodiments of the invention; 

Figure 50 shows steps in a method for processing 
data used for hydrocarbon .extraction from the earth, 
according to preferred embodiments of the invention; and 

3 0 Figure 51 shows further detail of steps in 

generating an efficient arid robust subdi vision of shapes, 
according to preferred embodiments of the invention. 
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DETAILED DESCRIPTION OF THE INVENTION ; 

While Low-level curvature-based methods as applied 
to implicit surfaces are relatively complicated to 
develop/ they do not suffer numerical stability problems. 
By comparison, on a smooth Riemannian manifold the 3rd 
derivatives of the square of the signed distance function 
describes the norm of the 2nd Fundamental Form and the 
mean curvature. For triangulated surfaces, this result 
is difficult to apply, because numerical evaluation of 
3rd derivatives is not guaranteed to be numerically 
stable. 

According to the invention, differential geometry 
methods of surface representation' will now be described:- . 
Many geoscience phenomena are related to some form of 
fluid flow. If the fluid phenomena under study involve 
surface tension, then mean curvature flow (mcf) and 
variants such as volume preserving mean curvature flow 
(vpmcf) are accurate modeling tools. (For those 
unfamiliar with mean curvature flow, imagine the fluid 
front moving normal to each, particle of the fluid front.) 
The behavior of mcf and vpmcf are well understood when 
either is applied to a smooth convex surface, star-shaped 
surface, a surface of rotation, or an entire graph. 
5 A reservoir structural framework does not seem to be 

an ideal input to mcf, because noise levels degrade the 
accuracy of analytical approximations and framework 
surfaces in general are not well approximated as convex, 
star-shaped, etc. This perception is erroneous. Our 
0 investigation of curvature-based modeling shows that mcf. 
can be used to semi-automate its own noise suppression. 
Given smooth surface data, according to the invention a 
method is provided to decompose that surface into a 
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connected sum of star-shaped or entire graph or 
axisyrametric surface patches. 

The mathematical foundation of mcf is substantial, 
so we seek a unified mathematical description of shape 
5 and its evolution. The concept of a fibre bundle is 

satisfactory. Nakahara presents a formal definition of a 
fibre bundle. Standard examples of fibre bundles are 
vector fields, e.g., velocity fields, and tensor fields, 
e.g., stress fields and elastic fields,' evaluated over a 
10 sub-volume. In classical differential geometry, curvature 
properties of surfaces are economically studied in a 
fibre bundle setting. We have found that a 4D fibre 
bundle representation of a reservoir framework is no more 
' difficult to write down than is a 3D fibre bundle " 
15 representations. The economy of mathematical 

representation is attractive now from the research view. 
It will fc>e attractive from the engineering view, since 
reuse of concepts limits the amount of technology that 
must be mastered. 
20' The fibre bundle representation of a surface of 

revolution has the following parts. 

1. A base that is homeomorphic to a loop, 
\ e.g., a circle, an ellipse, or a polygonal closed 
curve . 

25 2 . • A typical fibre that is. homeomorphic to a 

compact connected part of either a- conic section or 
a polygonal line. 

3 . A structure-group that is a group of 
dif feomorphisms that smoothly deforms any* instance 

30 of the typical fibre to any other instance. 
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Here are a few shapes that are frequently 
encountered in earth modeling represented as fibre 
bundles . 

a. A torus has a circle base and a circle 

5 . typical fibre. Its structure group is the rotation 

group SO (3) . 

b. A cyclide is a torus, but the structure 
group is the rotation group extended by the scale 
factor diagonal group. Intuitively, it is a skewed 

10 torus . 

c . Any compact planar region R with a 

polygonal boundary 3r has a natural bundle 
structure. To see this, compute the region's 
bounding box B and embed R in B. (Think of B as a 
15 solid.) If there exists a part of dR that does not 

intersect any of the box faces, then attach a skirt 
S of normal rays to 9. Then S, R, and 9b enclose a 
sub -volume V. 

We have found that V is the preferred bundle. The 
20 base of the bundle is 9r and a' typical fibre is a 

polygonal line. The region R is formed from the rotation 
of a line segment emanating from the centroid of R and 
joined to the boundary 3r. The length of the fibre 
changes instantaneously. The structure group is the 
25 Euclidean group of rigid body transformations. The 
following diagram summarizes this m construction. 

Figure 1 shows on the left an open surface, i.e. a 
.surface that does not -enclose space; on the right, a 
closed region (the large box) is subdivided by the same 
30 surface forming an internal boundary. Figure 2 shows a 
salt mass's steep flanks and overhangs, it also shows an 
example of the cyclide shape in depth imaging. This is an 
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example that challenges existing commercial software. 
WesternGeco commercial processing used to construct this 
velocity model has the common limitation of accepting 
- single z -valued ("height field") data only . According to 
5 the invention we describe below a set of planes in the 
volume of interest such that a multiple z-value body can 
be subdivided into sections such that each section is 
single valued with respect to one of the planes . In other 
words, each plane parameterizes a section of the, 

10 reference multi z-valued surface. Each of these planes 
are equipped with a rotation matrix and translation 
vector so that the application can orient the surface 
section so that it appears to be single valued with 
respect to one of the coordinate planes. (It may be that 

15 " the rotated and translated section is normal to the (x,y) 
coordinate plane, which is frequently unacceptable to 
grid-based applications . ) 

Material properties of a volume of earth can change 
during an evolutionary process. The time scale of the 

20 evolution can vary from wall clock to calendar to 
geological record. According to the invention, an 
evolutionary process is represented using a 
generalization of the fibre bundle method that is 
employed for .shape. This representation is called an 

25 "evolutionary process". Here are its components. 

1. The base is a line with period X B . A 
trajectory on the base serves as a clock. The 
trajectory's sampling increment has units that are 
30 appropriate to wall clock sampling, calendar 

sampling, or geological record sampling. 
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2. A trajectory has compact support. This 
means that the process is defined on a closed at the 
beginning and open at the end bounded interval of 
the helicoid. 

5 3 . Two evolutionary processes Pi and P 2 can be 

summed if and only if on a shared interval of the 
base helicoid Pi followed by P 2 is identical to P 2 
followed by Pi. The sum of Pi and P 2 is 
Pl(t) [P2(t) [v] ] where t is a point in time and v is 

10 a paint in the volume domain of interest, 

4. A fibre is a compact path- connected sub- 
volume. A fibre is associated to each base point. In 
other words the process time stamps every sub-volume 
that it affects. 

15 5 . A structure group is a 1 -parameter 

subgroup of dif f eomorphisms that define an 
evolution, , e.g. , vpmcf. Since the structure group 
consists of dif f eomorphisms the process must be 
invertible. In particular, the process cannot induce 

20 singularities in material space. It is preferable to 

support singularities, so we allow singularities to 
develop at the end of a process's time interval. 

We emphasize that we seek only a geometrical 
25 evolution - not necessarily the true physical evolution - 
that describes part of the structural framework. In the 
fqllowing example, we recognize uniformly expanding mean 
.curvature flow creates the shape of the individual 
layers. Each discontinuity in the sediment terminates the 
30 current flow model and is part of the initial conditions 
that define the new flow. The ensemble of flow problems 
describes the formation, but a more economical 
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representation can be defined if we can treat the entire 
set of restricted flow problems as a single mean 
curvature flow problem where evolution continues beyond 
the intermediate singularities. We lopk for an evolution 
5 of the. set of initial conditions for the individual flow 
problems/ given just the oldest sediment as an initial 
condition. 

We describe a sedimentary sequence as a 4D fibre 
10 bundle, 

1. The base<of the bundle is a fini.te length 
piecewise linear curve. 

2. A fibre is an infinitesimal layer of 
15 sediment. 

3. ' The bundle's structure group is a set of 
dif f eomorphism groups such that each group defines 
an instance of uniformly expanding mean curvature 
flow. Each leg of the bundle base defines a separate 

20 mean curvature flow problem. The transition between 

mcf problems is exactly matched by the discontinuity 
in the sediment . 

Figure 3 shows an example of a 4D representation of 
25 a field in Turkmanistan. Each layer is a distinct mean 
curvature flow where the discontinuity is a curvature 
flow terminator. Working backwards from a flow 
termination, we see that it is much easier to. identify 
the flow components in the image. Each flow component is 
30 a uniformly expanding set of solutions to- mean curvature 
flow (or volume-preserving mean curvature flow) . Ecker 
shows in his lecture notes that mean curvature flow can 
• evolve cracks and holes in a smooth background. See, K. 
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Ecker, Lectures on Regularity for Mean Curvature Flow , 
http : / /web . mathema tik . uni f reiburg . de/mi / analysi s / lehre/WS 
0001/Ecker_WS0001.ps. 

Therefore its usefulness is not limited to modeling 
5 smooth elastic behavior. 

According to the invention, we show how to construct 
a meandering river as a fibre bundle. Figure 4 is an 
image of the MacKenzie River Delta.. 

In Figure 4, the river does not intersect itself, 
10 i.e., no oxbow structures are evident. Also, the river 

appears to have constant width. The structure as a fibre 
bundle is clear. 

1. The bundle base is a line segment. 
15 2. The fibre is any convenient approximation 

of the river channel in cross section, e.g., a 
trapezoid! 

3 . The structure group is the extended 
Euclidean group, which consists of rigid body motion 
20 plus scaling. 

At each point along the channel we measure the cross 
section and record its position relative to the image 
coordinate axes whose origin is the lower left corner of 

25 the picture. Given the previous position of the cross 
section, we compute the update to the (x,y) plane 
rotation and translation. 

We turn now to the question of recognizing the 
elementary shape basis elements in an implicit surface 

30 representation. Recall that we represent an implicit 
surface as the zero level set of the signed distance 
function (sdf ) . Our constructor solves the sighed 
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distance function on a 3D structured grid with tri-linear 
interpolation as a local approximation to sdf . 

By definition tri-linear interpolation T(x,y,z) on a 
grid cell is 
5 , 

T(x, y,z) = A 0 + A 1 x + \y + A^z + A A xy + A 5 yz + A 6 xz + A 1 xyz* 

where the coefficients {A k } are defined on the grid 
cell corners. 

10 • Each term in T(x,y,z) is an independent 3D spherical 

harmonic polynomial in Cartesian coordinates. For 
convenience we enumerate these spherical harmonic 
polynomials, using the classical notation. (See R. 
Baerheim, Coordinate Free Representation of the 

15 Hierarchically Symmetric Tensor of Rank 4 in 

Determination of Symmetry , Ph.D. thesis, University of 
Utrecht, #159, 1998, Appendix, pg. 141-143) . 



Y(m,n) 


Monomial 


Y(0 f 0) 


1 


Y(0,1) 


z 


Real (Y(l; 1)) 


X 


Imag(Y(l,l) ) 


y 


Real (Y (2,1) ) 


xz 


Imag(Y(2,l)) 


yz 


Imag(Y(2,2) ) 


xy 


Imag(Y(3,2) ) 


xyz 



20 

Here are the associated symmetry transformations and 
the corresponding isomorphism groups involving these 
spherical harmonic polynomials. 



Term 


Symmetry generators 


Symmetry group 


Constant 


Constant 


SO(3) 


X 


[x] ~ [x] 


SO (2) 


Y 


[y] - [y] 


SO (2) 
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z 


[z] ~ [z] 


SO (2) 


X + Y. 


[x, y] ~ [y, x] 


Z/2Z 


X + x + Ai 


Y + Z 


[y, z] - [z, y] 


Z/2Z 


Y + Z + YZ 


X + 2 


[x, z] ~ [z, x] 


Z/2Z 


X + Z + XZ 


X + Y + Z 


[x, y, z] ~ [y, 

z, x] 
[x, y] ~ [y, x] 1 


Tetrahedron 
group of order 
12 


XYZ 


[x, y, z] - [y, 
z, x] 
[y, z] ~ [-y, -z] 
[x, y] - [y, x] 


Symmetric (4) 


XY 


[x, y] ~ [-x, -y] 


Z/2Z 0 Z/2Z 


[x, y] - [y, x] 


YZ 


[y, z] - [-y, -z] 


[y, z] - [z, y] 


XZ 


[x, z] - [-x -z] 


[x, z] - [z, x] 


XY + XZ 


[x, y, z] ~ [-x, 

-y / -z] 
[y, z] - [z, y] 


Z/2Z 0 Z/2Z 


XZ + YZ 


[x, y, z] - [-x, 

-Y/ -z] 
[x, y] - [y, x] 


XY- + YZ 


[x, y, z] -[-x, 

-y, -z] 
[x, z] ~ [z, x] 


XY + YZ + ZX 


[x, y, z] - [-x, 

-y, -z.] 
[x, y, z] - [y, 

z, x] 
[x, y] - [y, x].. 


Symmetric (4) 



Symmetry analysis is a form of spectral analysis 
applied to the discrete spectrum that is associated with 
spherical harmonic polynomial expansions. Our estimates 
5 of tri-linear interpolation coefficients are noisy, so we 
need a threshold for dismissing spectral lines. 

Figure 5 shows some combinations of involved 
spherical harmonic polynomials, presented in spherical 
polar coordinates. 
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We are interested in conic section fibre bundle 
shapes. Here is the correspondence of shape to symmetiry 
group . 



Symmetry group 


Candidate shapes 


SO(3) 


Sphere 


0(2) 


Cylinder, spherical torus, 
elliptical torus, hyperboloid 


i SO(2) 


Cone, paraboloid 


Symmetric (4) 


Cube 


Dihedral (4) 


Ellipsoid 


i Dihedral ( 3 ) 


Triangular prism 


z, ©z, 


Rectangular prism 




Cyclide 


{1} 


Noisy data 



* According to the invention, mean curvature flow and 
volume -pre serving mean curvature flow will now be 
discussed in. further detail. We define a few. terms that 
appear frequently herein. Given a 2D manifold M and a 
10 point p e M, let {pi, P2> be a local coordinate system for 
a region containing p and let n be an outward pointing 
normal at p. Finally the Euclidean inner product of 
vector Ui and Vj is denoted by <Ui, Vj>. 

Definitions 

The . Fir st Fundamental Form (1 st FF) is g u = 

15 The inverse of the 1 st FF is denoted by g iJ . 

The Second Fundamental Form (2 n ? FF) is A = [h ( j] = " ^2 p y" >n )' 

The Weingarten map is W/ = g ik h H . {Einstein notation used.).. 

The eigenvalues of the Weingarten map are the 
principal curvatures . * The trace of the map ■ is the mean 
curvature , the determinant of the map is the Gaussian 
20 curvature . The norm of the 2 nd FF |a| 2 is defined as the 
sum of the squares of the principal curvatures. 
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We define mean curvature flow and volume-preserving 
mean curvature flow. 

Notation 
5 . 

(MJ is a family of evolving smooth 2D manifolds such that 

M = M„ is given. 

N(x,t) is the normal at x e M t . 

H(x,t) is the mean curvature evaluated at x 8 M t . 
10 h(t) is the average value of H(x,t) on M t . 

Mean curvature flow (mcf ) is defined as the solution 
to the initial and boundary value problem 

<Ms- = -H(x,t) -N(x,t). [4] 
dt 

15 Volume -pre serving mean curvature .flow (vpmcf ) is 

defined as the solution to the initial and boundary value 

problem f^ = -(ff(jc,0-/i(0) -N(x,t). [5] 

Figure 6 illustrates that, except in ideal 
20 circumstances (when the input is a smooth closed convex 
region or a graph) , a surface evolves under mcf to a 
point. Technically, mcf develops a singularity in finite 
time. See, K. Museth, D. Breen, R. Whi taker, A. Barr, 
Level Set Surface Editing Operators , SIGGRAPH 2002. 
25 A straightforward way to prevent annihilation of 

this shape is to stop the mcf after some number of time 
steps, inspect the results, and maybe resume the process. 
This is not convenient in a production environment. In 
fact,, the" manner in which uncontrolled mcf annihilates a 
30 shape enables us to attach a recognition procedure that 
mcf can call to decide when to ask the human for 
permission to resume the figure's evolution. 
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Figures 7a- f forms a series of six images showing 
vpmcf suppressing noise. See, A. Kuprat, A. Khamayseh, D. 
George, L. Larkey, Volume Conserving for Piecewise Linear 
Curves, Surfaces, and Triple Lines, Journal of 
5 Computational Physics, 172 (2001), pg. 98-118.. In Figure 
7a, the southern hemisphere is corrupted with a 
significant amount of Gaussian noise. By Figure 7c, it 
makes sense to consider how much additional noise, if 
any, must be removed before the smoothed result is an 

10 acceptable approximation. 

( We are interested in surfaces of revolution. The 
following result has been obtained regarding the 
behaviour of surfaces of rotation under vpmcf. See, M. 
Athanassenas, Volume-preserving mean curvature flow of 

15 rotationally symmetric surfaces . Comment. Math. Helv. 
72(1997), pg. 52-66. 

. Theorem (Athanassenas) 
20 Let Mo be a smooth rotationally symmetric surface 

enclosing a sub-volume V. Let S be a slab of thickness x 
> 0 that is bounded by the z = 0 and z = x height field 
planes. Suppose that M 0 satisfies the following two 
conditions . 

25 

a. Mj^has a line of intersection in the z = 0 and z = 
x height field planes. 

b. The end points of each line of intersection 
between M 0 and S meet an end of Mo at a right- 

30 angle. 

If the volume of V is greater than x times the area 
of M, then M t evolves under vpmcf to a cylinder. 
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As an example of the condition on the lower bound on 
volume, consider the case of a right cylinder of height T 
and radius r. This cylinder has volume V equal to 7tr 2 T, 
while its surface area M is equal to 2nri. Hence the lower 
5 bound on the volume is satisfied exactly when 7ir 2 x > 

2nxi?,x.e, r > 2x. Heuristically , the bound is satisfied 
for w squat" cylinders. We observe that a sphere cannot 
satisfy the volume to area relationship. To do so would 

imply that there exists r > 0 such that ' —^r 3 /4^r 2 =->2r . 

I'O Figure 8 illustrates the orthogonality condition of 

the theorem proposed by Athanassenas . In this diagram we 
show two vertical cross-sections. In the figure on the 
left the intersection of the figure .with the top and 
bottom planes must satisfy the right angle hypothesis. 

15 The curvature flow takes care qf the irregular vertical 
surfaces either end of the figure. In finite time vpmcf 
generates from the figure on the left the figure on. the 
right. 

Volume preservation is essential for earth model 
. 20 applications, so we prefer vpmcf to ordinary mcf . Later 
in the description/ we will use this theorem to reduce 
the discrepancy between an idealized representation ■ of 
shape and a sampled data surface. 

Another important, class of smooth surfaces are those 
2.5 that are star-shaped. A surface is star-shaped if there 
exists a point P on the surface such that a line segment 
EQ. that is entirely, contained in the surface can join 
.every other point Q in the surface. 

It has been shown that star- shaped closed smooth 
30 surfaces are stable under mcf. See, K. Smocyzk, Star- 
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25 



shapfed hypersurf aces and the mean curvature flow , 
Manuscripta Math. 95 (1998), pg. 225-236. 

We show now thatvpmcf suppresses additive high 
frequency harmonics before decaying the underlying low 
frequency shape signal. Vpmcf enjoys the property that 
surface area is always decreasing in time. 



£^=- L (*-*)^,,o 



[6] 

dt 

10 We frequently need to compute mean and Gaussian 

curvature for a single valued surface over a plane, i.e., 
a graph. Suppose that S = S(x,y). Then . the formulae for 
mean curvature H and Gaussian curvature K are as follows. 



15 ^ = "i TZT and M ~ 7T~ 1 , \3 

a+s x 2 +s v 2 f 2 i.KZT^ 



We model a surface as a collection of patches, where 
each patch is taken from a surface of revolution, say S = 
(<|>(v) *cosu, <|>(v)*sinu, y(v)), such that v e^a/b), u 
20 6 (0,271)/ and \|f(v) * 0. The formulae for mean curvature. H 
and Gaussian curvature K for- a surface of revolution is 
as follows 



L // 



d 2 0 



For the proofs, see, M. do Carmo, Differential 
Geometiry of Curves and Surfaces, Prentice Hall, 1976, 
pages 1 '62 -163 . 

We need to understand how mcf behaves when applied 
3 0 to a noisy surface S, where S is either single valued 
over a plane or is a surface of revolution. The above 
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formulae imply that if S is an entire graph, then the 
harmonic number N scales the mean, curvature. Similarly, 
if S is a surface of revolution, then N 2 scales the mean 
curvature. (Spherical harmonic polynomials have a cosine 
5 factor, so these estimates remain valid under the natural 
heat equation ei gen- function expansion.) Consequently, 
the mean curvature integral involves a scale factor of 
either N 2 or N 4 . Therefore the rate of change of surface 
area becomes very negative as N increases . We conclude 

10 that vpmcf eliminates high frequency harmonics during the 
initial stages of the evolution ahead of low frequency 
harmonics that form the ideal shape. A maximum or minimum 
in the ideal .shape will appear to move as vpmcf 
eliminates corrupting noise, but the critical point's 

15 location will stabilize as the noise disappears. 

Importantly, when a critical point's location stabilizes, 
it does not move later in the evolution unless vpmcf 
begins to decay the shape. We remark that mcf can be 
turned off locally by resetting the mean curvature to 

20 zero. 

We turn now to the question of how to decide when 
noise suppression turns intq shape decay. Refer again to 
Figures 7a- f, which show the smoothing of a noisy sphere 
under vpmcf. Athanassenas ' s theorem does not apply, since 
25 the sphere fails the volume bound assumption. However, it 
still makes sense to apply Kuprat et al.'s vpmcf 
procedure to the shape. When we do this we get Figures 
7c-f. 

It is reasonable to say that noise suppression 
30 reduces Figures 7a-b to Figures 7c-d. It is harder to say 
if image noise suppression or shape decay accounts for 
the transformation to Figures 7e-f . Therefore, we want 
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vpmcf to proceed without manual intervention to eliminate 
noise but to seek guidance when the flow causes the 
under lying shape to decay. Here is a way to monitor the 
elimination of noise and detect decay of the critical 
5 points in an underlying shape that relies on vpmcf. 

1. Given a characteristic length, we embed a series of 
conical shapes in the surface. The base of the cone 
is either the outermost ID boundary of the surface 

10 or it is an inflection point curve surrounding the 

patch that contains the maximum or minimum. 

2. We recognize minima and maxima in its 2D enclosure 
using a standard nested bounding box algorithm. If 
the surface is implicit, then we can localize the 

15 region by testing the implicit surface's 

• triangulated parent. 

3 . Use this 2D enclosure as the top of the (noisy) 
cone. Note that we cannot assume that each maximum 
or minimum is topologically isolated, e.g. 1 , a ridge 

20 might exist. 

4. Apply any mean curvature flow procedure, e.g.,' D3b, 
to the surface and monitor the effect on the cone's 
slant height. Mean curvature flow theory says that 

. the slant height decreases as the noise is removed 
25 (it becomes straighter) . However, when the mqan 

curvature flow begins to. destroy the intrinsic shape 
of the surface in the vicinity of the cone, then the 
slant height will decrease beyond a threshold. 
Therefore, the test for adequate smoothing is to 
30 compare the ratio of slant height before smoothing 

to the result obtained after each iteration. 

Figure 9 is an aerial image of part of Big Bend 
National Park, showing the approximation of a. plateau to 
35 a characteristic length scale cone. It is difficult to 
precisely locate the maxima on the plateau, but it is 
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4 

easy to enclose the region containing the maxima in a 
tight loop. It is unimportant that the cone does not have 
a circular cross section. 

Figure 10 is a side view of a noise cone structure. 
5 When monitoring noise suppression on the surface of a 
geobbdy, e.g., it is usef ul . to have a noise-monitoring 
device for thin undulating cross section. We replace a 
torus by a cone. We use the. cone to monitor shear 
stretching and erosion of an interface. A clear instance 
.10 of this phenomenon is Figure 7c. 

According to the invention, mean curvature flow and 
framework I/O will now be discussed in further detail. 
Singularities mark an end to the smooth evolution of a 
15 shape under mcf . A singularity is frequently easier to 
identify on an image than is the interface that 
corresponds to the precise start of the flow. The flow 
imposes a natural partition of the framewQrk. Each region 
in the partition is a* self-contained expression of mcf. 
20 When we think about sending and receiving an update to a 
framework, we prefer to send and receive a mcf problem 
with boundary conditions rather than an opaque byte 
stream. We specify the solution form - whether the flow 
uniformly expands or contracts the solution or smoothes 
25 the initial data -and provide beginning and end surface 
data. We describe the intermediate surfaces by recording 
the time states that correspond to the intermediate 
surfaces . We also send critical point data for each 
intermediate surface. Interested applications run the 
3 0 vpmcf script and reconstruct the. update locally. We trade 
fast CPU for less fast. I/O. 



24 



WO 2004/057375 




PCT/GB2003/005395 



Here are three illustrations of this idea. Figure 11 
is a NASA satellite 'image of the Labrador Trough. In this 
image we notice that the sediment resembles a 
longitudinal cross section of a brain with an attached 
spinal cord. Singularities separate "brain tissue" from 
the stem of the w spinal cord". A singularity also 
subdivides the * cerebellum" . 

Figure 12 shows s a sequence of folded sediment on 
the coast of the Gulf of Oman. We model this as a 
uniformly expanding (seen right to left) solution to mcf . 
Ecker proved that a uniformly expanding solution to mcf 
is given by the equation M , = • M x , provided that the 
initial fold is approximately a cone. 

Figure 13 is an image illustrating progressive 
flattening of an overburden covering a large salt 
intrusion. We notice that the fault block that occupies 
the left half of the image (indicated N by the lower white 
arrow) displays a gradual flattening in space and time as 
the top of the salt's intrusion is progressively reduced 
as the sediment was deposited on top of the salt 
(indicated by the upper white arrow) . 

According to the invention, a preferred technique of 
2D parameterization will now be discussed in further 
detail. Many surface editing operations are more 
efficient when the operator can access the surface's 
parameterisation. Height field data uses the (x,y) co- 
ordinate plane as the parameter space and coordinate 
projection as the parameter space mapping. We explain how 
to obtain a parameterisation of a smooth surface S with 
no assumptions regarding the orientation of the surface 
relative to the (x,y) .co-ordinate plane. 

25 
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We claim that there exists an invertible projection 
of S onto a cylinder that has 0, 1, or 2 end caps. 

1. Given any two closed curves C A and C 2 on S. 
We partially order the two curves, say C x < C 2 , if 

5 there exists a 2D region R, on S such that dl^ = C 2 

and C x C R, We call the 2D region of S that is 
contained between *C 2 \ C x the collar that is bounded 
by C 2 and C 1# 

2. From the classification of smooth surfaces 
10 in R 3 we know that S has Euler Characteristic %(S) 

equal to 0, 1, or 2. The Euler Poincare Theorem says 
that x(S) = # (maxima) + # (minima) - # (saddle 
points) . Note that we use height as our Morse 
function, so the term "critical point" is a 
15 contraction for "height field critical points 

3 . We will prove our claim by induction on 
the number of saddl'e points on S . 

4. Suppose that the number of saddl'e points 
on S is zero, 

20 a. If %(S) = 2, then Reeb's Theorem says 

that S is dif f eomorphic to a round sphere, 
which is dif f eomorphic to a cylinder with two 
end caps . 

b. If JC(S) = 0, then S must be 

25 dif feomorphic to a plane which is dif f eomorphic 

to a cylinder with one end cap. 

c. If x(S) = 1* then S has either 1 
maximum or 1 minimum but not both. 

d. Construct a 1 st order Taylor series 
30 with. center equal to the critical point and 

quadratic -order error term. 

e. Notice that the gradient term is 
trivial, since we expand about a minimum or 
maximum. 
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f . We apply a translation to the surface 
so that expansion point is the origin and has 
height equal to 0. 

g. . We conclude that after a suitable 

5 ■ translation and (x,y) plane rotation a 2 term 

Taylor expansion about the critical- point with 
a quadratic error estimator simplifies to a 
pure quadratic term. 

h. Given 8 > 0, we fit the Taylor 

10 expansion such that the error estimate is less 

than £. 

i. Then the Taylor series expansion 
within C 0 is a cap /cup. We know that a cap /cup 
is dif f eomorphic to a plane and that a plane is 

15 dif f eomorphic to cylinder with one end. cap. 

Furthermore; the open end of the cylinder has 
boundary C 0 . . 

5 . Suppose that the number of saddle points 
20 on S equals 1. 

a. If %(S) = 2, then S possesses 
either. (i) 2 maxima and 1 minimum or (id) 2 . 
minima and 1 maximum or (iii) 3 maxima or (iv) 
3 minima. 

25 b. Case (iii) is impossible, because 

Morse Theory says that a maximum contributes a 
2D cell to the shape. If case (iii) were true, 
then 

c. Case (iv) is impossible, because 
30 Morse Theory says that a minimum contributes a 

0D cell to the shape. If case (iv) were true, 
then S contains ho "2D cells. 

Figure 14 is a diagram illustrating the Morse 
35 theoretical cell decomposition for a simple configuration 
of a capped and bent cylinder. 
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Morse Theory explains the cell decomposition of the 
shape. (See M. Hirsch, Differential Topology, Springer 
Verlag, 1976, pg. 156-164.) 

• d. This shape has a saddle point at 
location C. 

i. Morse Theory says that the 
' surface contains a ID cell that is 

attached to location C. 

ii. We observe that this ID cell is 
not a loop. Specifically, the point at 
location D is not part of the ID cell. 

e. This shape has a minimum at location 

D. 

i. Morse Theory says that the shape 

contains a 0D cell that is attached to 
location D. 

f . This shape contains two maxima at 
locations A and B. 

i. Morse Theory says that the shape 

contains two 2D cells, and that each 2D 
cell is attached to a maximum. 

ii. The boundary between the two 2D 
cells is the union of the 0D cell that is 
attached to location D and the ID cell 
that is attached to location C. 

iii. We decompose each 2D cell into 
the union of a Taylor series cap expanded 
about the maximum and a collar that is 
bounded by^ the loop suspended between 
locations C and D and the -Taylor series 
circle of convergence. 

g. Suppose that %(S) = 1. The Euler- 
Poincare formula says that the surface' contains 
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either (i) 1 maximum and 1 minimum or (ii) 2 
maxima or (iii) 2 minima. 

i . An example configuration for 

case (i) is the upper half of a torus. 
5 ii . An example configuration for 

case (ii) is a flat surface with two 
hills. 

iii. An example configuration for 

case (iii) is a flat surface with two 
10 valleys . 

h. Finally suppose that %(S) = 0. Then S 
contains either 1 minimum or 1 maximum, but not 
both. An example of this configuration is a 
15 single mountain or .sinkhole attached to a 

plane. 

6. We have specified the parametric mapping 
when S contains either 0 or 1 saddle, point. 

20 

7. Now suppose that contains N-l saddle 
points." Choose any saddle point and pass a height 
field plane through it, subdividing S in two. At 
least one of the two parts contains fewer saddle 

25 points than does S. Define S* to be this part of S. 

a. By induction we know that S* is 
dif f eomorphic to either a cylinder with a hole 
in its. lateral surface, or a cylinder that has 
0 or 1 end caps. 

30 b. Closing the hole in S* corresponds to 

either closing a hole in the walls of a - 
cylinder or to attaching an end cap to the 
cylindrical surface that we created in .xyz 
space . 
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Figures 15-17 show diagrams to aid in the 
under standing of the bulls eye construction. In Figure 



29 



WO 2004/057375 




PCT/GB2003/005395 



16, we invertibly map the possibly patched A to the 
planar unfolding of its interior facetted cylinder. 
Figures 17a and 17b illustrates how to perform the bulls 
eye construction in the region of a saddle point. The 
5 saddle point is indicated with the dot 10 in Figures 17a 
and 17b This process is preferably repeated for every 
saddle point of the shape. 

According to the invention, tri-linear 

10 interpolation, implicit surfaces, and critical points 

will now be described in further detail. In particular 
we herein discuss the analysis of the tri-linear 
interpolation approximation to the signed distance 
function on a 3D structured grid. We derive necessary two 

15 conditions that the tri-linear interpolator must satisfy 
if the grid cell contains a height field critical point. 
One condition specifies a relationship between the tri- 
linear coefficients and the critical value. The other 
condition establishes a second relationship among four of 

20 the tri-linear coefficients. Taken together, we show that 
a height field critical point for an implicit surface 
must be embedded in the grid cell faces and can never be 
found in the interior of the grid cell unless the 
implicit surface is a plane. 

25 Let G denote a 3D structured grid with typical grid 

cell C. We assume that an implicit surface S is contained 
in G and in particular that C nS 0 . We denote the tri- 
linear interpolation function on C by a and we assume 
that C is small enough that o is a good approximation to 

30 the signed distance function that implicitly defined S. 
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We begin with' the definition of tri-linear 
interpolation 

• <T(x,y,z) = A(x,y)+ B(x,y) z (1) 

5 

where 

A(x, y) = a 0 + a x x + a 2 y + a 3 xy 

B (x,y) = b 0 + b x x + b 2 y + & 3 jc)>. {2) 

10 For later reference, we note that A and B are 

harmonic , i.e., 
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V 2 A = V 2 £ = 0 < 3 > 

We are interested in the zeroth level of a. If a = 

0, then either B =0 or B * 0- We suppose that B * 0. 
Then 

20 Z D/ \' (4) 

We know that height field critical points are found 

when 
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dz _ dz 
dx dy 



25 — = ^-=0. ( 5) 



This condition is satisfied when 



B A ^B = o 

ox d x 

(6) 

=0- 

o y d y 

Rewriting (6) 
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8 A ^ A 



A 




8 X 




8 


y 


B 




3 B 




3 


B 



(7) 



a y 



Using (7) 



3 y 9x 3 x 

8 A 95 3A . . _ 

(— .-= -r ) + ( B a 3 - A b 3 ) - (8) 

By ox ox o y 

( B a 3 - Ab 3 ) = 0. 

The Hessian matrix of 2 nd partial derivatives 
evaluated at the critical point is 



3 2 Z ^ " dx K " 3* " 8.x J _ 



a* 2 a 4 



= o 



3 * _ 3y 9y = Q (9) 

dy 2 5 4 

a 2 A , a 2 s 



d 2 z a 8;c8y 9jc8y = Ba 2 -Ab 3 _ Q 
dxdy B 4 ~ ~ B 4 

We conclude that the Hessian matrix of 2 nd partial 
derivatives is identically zero. We are led to consider 
what kinds of shape have the property that in a local 
neighborhood of a critical, point the matrix of second 
partial derivatives is identically . zero*. 

We have enough information to compute the Gaussian 
curvature of z = z(x,y) . Recall that for a graph z = 
z(x,y) that its Gaussian curvature K(x,y) is 
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B 2 z B 2 z , B 



z 



10 



20 



30 



v r v „\ d* 2 dy 2 BxBy 

K(x,y) = (10) 

dx By 

We notice that the signum of K(x,y) depends only on 
its numerator. Substituting' (-3) into the numerator of 
(9), we discover that the numerator of K(x,y) is 



< dz 2 ^2 _ (n n A h ,,M dB BA BB. 
W -~ iB ' a ^ A ^ ) + ( By-'^ + ^^y- 



Hence the Gaussian curvature K(x,y) < 0 everywhere. 
As an illustration we consider the "monkey saddle " 
M(x,y), which coincidentally is also a spherical harmonic 
polynomial M(x,y) = Re(Y 33 ) = .x 3 - 3xy 2 . The monkey saddle 
has a single critical point, which is the origin (0,0). 
15 The matrix D of '2 nd partial derivatives is 



(12) 



6x ~6y 
-6y -6x 

Substituting (x,y) = (0,0), we find that the 2 nd 
partial derivatives matrix is identically 0. We conclude 
that the origin is a degenerate saddle point for M(x,y) . 

Figure 18 shows two views of an example of a monkey 
saddle. In particular the monkey saddle is M(x,y) given 
by the equation above. The image on the left is shaded 
25 according to mean curvature, while the image on the right 
is shaded according to Gaussian curvature. We have placed 
a white dotted circle around -the critical point.- - - ' 

The critical point is a degenerate saddle point. We 
see this by substituting the samples defined by the white 
circle and observing the pattern of +/- signs. As an 
example, let 0 < a < b. Then ((+/-) a, ( + /-) b) produces the 
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pattern (-,-), (-,+), (-,-)/ and (-,+) clockwise from 
Quadrant I. 

The Gaussian curvature of the monkey saddle is 
• (Gray, pages 382-383) 

5 



^/ 36(x 2 + y 2 ) . • (13) 

K{X * 30 " ~(l + 9* 4 +18*V + 9y 4 ) 2 

(See, A. Gray, Modern Differential Geometry of 
Curves and Surfaces with Mathematica 2nd Edition, CRC 
10 Press, 1998, pages 382-383). 

Summarizing the analysis for B ^ d, we have shown, 
the following. 

1. A grid cell contains at most .one .height field 
15 critical point. 

2. If a grid cell does contain a critical point, then 
the critical point's coordinates are uniquely 
determined from the tri- linear interpolator's 
coefficients. We use this as a quick test for 

20 assessing the existence of a critical point in a 

grid cell. 

3. The matrix of 2 nd partial derivatives is 
identically zero when evaluated at the critical 
point . 

25 4. The Gaussian curvature K(x,y) < 0 and equals zero 

exactly at the critical point. 
5*. z = z('x,y) is harmonic, so by the Maximum 

Principle for harmonic functions on compact sets, 

we know that the minimum and maximum value for z 
30 on a grid cell will be found on the grid cell 

corners . 



34 



WO 2004/057375 



PCT/GB2003/005395 



It has been remarked that the behavior of the tri- 
linear interpolation function is determined in the region 

^ dz dz . dz 

5 of a critical point by the behavxor of — 9 —,and — . See, 

G. Weber, G- Scheuermann, H. Hagen, B. Hamann, Exploring 
Scalar Fields Using Critical Isovalues, http: 
/ /graphics . cs.ucdavis . edu/~hamann/WeberScheuermannHagenHa 
mann2002.pdf ("Weber"). We agree,- which is Morse Theory 

10 should be applied. 

Weber shows that a tri-linear interpolation function 
can have a critical point on the edge of a grid cell. if 
and only if the tri-linear interpolation function is 
constant along the edge. Weber does not assume that the 

15 tri-linear interpolation approximates a signed distance 
function. Weber also obtains a simple test for the 
existence of a maximum at grid cell vertex (0,0,0) by 
looking at the value of the function at the tetrahedral 
comers (1,0,0), (0,1,0), and (0,0,1). 

20 In our situation Weber's condition says that a 0 > max 

(ai, a 2 , b 0 ) implies that (0,0,0) is a maximum of T(x,y,z) 
= A(x,y) + B(x,y)*z. We note that T(0,0,0) = a 0 and that 



dT 



0 



dx 



(0.0,0) 




S(0,0) = fr 0 



0 



(14) 



dz (0,0.0) 



25 



Checking 2 nd mixed partials 
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10 



20 



dx 2 dx 2 dx 2 

dy 2 By 2 Z 'dy 2 

^ = 0 

dz 2 



3 2 T 


3 2 A 


3 x 3 y 


3 x 3 y 


3 2 T 


3 B 


d x d z 


. 3 x 


3 2 T 


• 3 B 


3 y 3 z 





B 



= a 



3 



d x d y 

ft ! + ^ 3 ? (15) 
2 + ft 3 x 



Therefore the first order Taylor series expansion of 
the tri-linear interpolation function about (0,0,0) 
inside a grid cell cube is given by the expression 

v 

T U, y, z ) = a Q + b 0 z + (b x + b 3 y)xz + (ft 2 + b z x)yz 

= a 0 + b 0 z + b l xz + b 2 yz + b 3 xyz + b 3 xyz (16) 

= a 0 + B(x, y)z + b 3 xyz 

This says that in a local neighborhood of (0,0,0) 
that the tri-linear interpolation function is a harmonic 
polynomial and therefore the Maximum Principle applies. 
15 We conclude that on the tetrahedron formed by the four 

grid cell vertices that maximum occurs at the grid cell . 
vertex a 0 or b 0 . We invoke Weber's assumption and conclude 
that the maximum occurs at a 0 . It is not necessary to 
assume the other inequalities in Weber's assumption. 



We consider now the case that a = 0 and B = 0. We 
begin by observing that o = 0 implies that B = 0 if and 
only if A = 0. Therefore, 
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y = . a o + a i x = b 0 +b J x 

a 2 + a 3 x • £ 2 +£ 3 jc* (17) 

Again we look for solutions to the zero gradient 
equation. 

dy a { a^ - a 6 a, ^ 

-7— = - - — ~r* = 0 implies that a,a y - a 0 a q = 0 

<ix (a 2 + a 3 jc) 2 • 1203 

. and {18 j 

^ " ~~ lb*+b*x) 2 = ° //n ^ Zi " that b x b 2 - fr 0 6 3 = 0 

Again we obtain a relationship among the tri-linear 
interpolator coefficients that must be satisfied in order 
that a grid cell face contain a critical point on a curve 
running through the grid cell- face. 

Suppose that a grid cell face contains a critical 

point of a curve. Again the second derivative — ~ = 0, so 

dx 2 

15 we characterize the nature of the critical point by 
sampling the gradient on a tight neighborhood of the 
critical point. 

We conclude that there exists at most 1 critical 
2 0 point on a curve running through a grid cell face and 
that the (x, y) coordinates of a critical point are 
uniquely determined in terms of the tri-linear 
interpolator's coefficients. 

25 According to the invention shape index and shape 

identification techniques will now be described in 
further detail. According to a preferred embodiment, 
] shape identification depends on a dimension - independent 
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measure of the principal curvatures at a point known as 
the Koenderink and van Doom shape index si. 



2 K . + K o 
s i = — ; a r c t a n — — = 

71 K 1 - K 2 



2 H 
— arctan — , = , 

x ^JH 2 - K 



10 



where fC L >fC 2 . 

Shown below is a chart of the shape index map's 
range [-1.0, +1.0], which goes from most concave to most 
convex. We' note that a shape index value of zero 
corresponds to a zero mean curvature surface, which for 
the case of a compact surface, equates to a catenoid or a 
compact planar region. 



Local shape 


si interval 


spherical cup 


1-V7/8) 


trough 


[-7/8,-5/8) 


rut 


1-5/8,-3/8) 


saddle rut 


[-3/8,-1/8) 


saddle 


[-1/8,1/8) 



Local shape 


si interval 


saddle ridge 


[1/8,3/8) 


ridge 


[3/8,5/8) 


dome 


[5/8,7/8) 


spherical cap 


[7/8,1] 



Cantzler et al. have computed a correspondence 
between shape attributes defined by a Gaussian and mean 
15 curvature value pair and the shape index map range. See, 
H. Cantzler/ R. Fisher, Comparison of HK and SC curvature 
description methods , 

http: / /www. dai .ed.ac.uk/homes/rbf/hc3dim.ps.gz. See, the 
following tables. 
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K<0 


K=0 


K>0 


H<0 


Saddle Wley 


Concave 


Concave 




• (SvHy) 


Cylinder (-Cy) 


Ellipsoid (-E1) 


H = 0 


Minimal 


Plane 


Impossible 




(MHy) 


(PD- 




H>0 


Saddle Ridec 


Convex 


Convex 




(SrHy) 


Cylinder (+Cy) 


Ellipsoid (+EI) 



Shape 


Index range 


Concave Ellipsoid (-E1) 


S€ [-V5/8) 


Concave Cylinder (-Cy) 


S € [-5/8,-3/8) 


Hyperboloid (Hy) 


S € [-3/8,3/8) 


Convex Cylinder (+Cy) 


S€ [3/8,5/8) 


Convex Ellipsoid (+E1) 


S € [5/8,1] 



As an example, we compute the shape index of a 



torus 



5 



10 



c . 



15 



20 



The torus 's parameterization is 

((a + fc-cosv)cosw, (a+£>-cosv)sinw, b-sinu), 
where we assume that a > b > 0. 

cosv 



a + b cosv 



The principal curvatures are 
Therefore - 1 — — = -(1+ — cosv). 

h-K a 

This quotient covers the range [-3.0, +1.0] 
the shape index covers the interval [-0.80, 
which corresponds to "trough" to "ridge" . 



and k 2 =— — . 



hence 
+0.50] , 



25 



Certain' principal curvature pairs are distinguished. 

1. If = K x , then the shape index is +/- 1.0, as 
the signumdCj) = -/+1. 

2. If 0 = Kj > iq, then the shape index is 

3. If Kj > iq= 0, then the shape index is 

4. If = -lq, then the shape index is 0. 

5. The shape index of a maximum is 1.0. 

6. The shape index of a minimum is -1.0. 

7. If the shape index is discretized in increments 
of size As , then the principal curvatures will 
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be discretized according to the following 
formulae 

Shape index expresses a relationship between 
principal curvatures. For example, we expect the shape 
index in the region of an inflection point is 
approximately 1/8. Plugging into the shape index 
definition, we discover that 



10 /r 2 = K L 



( , ,n 1 , , ^ 
tan ( ) - 1 

2 8 ' 

tan( — ) + 1 
V 2 8 



= -0.99023 K t 



The shape index enables a human to express a 
threshold change in curvature in a dimension- independent 
manner. This correspondence is fundamental to the 
15 robustness of our shape identification scheme. 

Measurements are always tainted with noise. Therefore it 
preferable to identify intervals rather than point values 
for attribute, correspondence. 

Definitions 

Let A,B £S and let s A ,s B be the shape index at A,B. 

20 By definition ten(—-s A ) = 0 A and tan(— • s B ) = 6* , where ' 

2 2 

A \ ty A J 



6 = where K A and K 2 A are principal curvatures at A. 

Likewise for 0 B . 



TIius tan(~(j B -s A y) = ?* B 0 * or (s B -s A ) =-arctan(-^ 



2 v \" 1-0*. 0* ~ v " " ' n ™\-0*.0 A 
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This formula relates the shape index principal 
curvature quotient in shape index increments. For example 

if {s B s A )=\, then tan(^) = 0.0034= ^^ . 

Holding 6 A fixed, then tan(^) • (1 - 6*6 A ) = 6 B -6 "implies that 

16 

tan(jjr) + 0" 0.0034 + g* 

= 1+tan( 2L).^ = 1 + 0.0034 <>>■ 
16 

Now we can utilize the formula relating the value of 
a shape index to the ratio of the constituent principal 
curvature . 

According to a preferred embodiment of the 
invention, techniques for curvature- adaptive sampling of 
a smooth surface will now be described in further detail . 

15. Due to popularity of 2D FFT methods, it is common 

practice to sample a surface based on a fixed spatial 
step size. For smooth data, we know that a low order 
Taylor polynomial expansion can supply sample data to 
whatever a priori precision and at whatever spacing is 

20 acceptable. There is no need to store a large array when 
an algebraic expression can be evaluated on demand. Our 
' method for curvature-based sampling is straightforward. 

1. We are given a cloud of points. Using Hoppe's 
25 tangent plane fitting procedure (Hoppe, pages 

22-25), we loop through the input point cloud 
fitting tangent planes. 
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2. At each tangent plane, we determine if the 

parent point is a height field critical point. 
We record the set of all critical points {C k } . 
We use the symbol T(C k ) to denote the tangent 
plane at C k 

3 . * At each critical point C k , we define a 

neighborhood N(C k ) about C k in the point cloud 
such that the tangent plane T(C k ) is a 
parameterization of the triangulation of the 
points in N(C k ) . Abusing notation slightly, we 
use N(C k ) to denote the triangulation of the 
point cloud. 

4. Using equations (4.7^ and (4.7 2 ) we estimate 
both principal curvatures and the associated 
shape index on N(C k ) . 

5. Given a shape index increment, we apply the 
analysis developed in Section 8 to partition 
N(C k ) into shape index equivalence classes. We 
chose to seed the shape index equivalence 
classes with critical points so that there is 

. the least ambiguity in the shape index 
assignment. 

6. If the critical point neighborhoods do not 
cover completely the implicit surface, then 
search the remaining surface for regular points 
that have an unambiguous shape index interval 
assignment . 

7 . Bin the remaining unassigned sample points to 
new shape index equivalence classes . 

8 . Assign each shape index equivalence class to a 
distinct grid cell. 
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9 . Done.. 

This algorithm uses just curvature information, 
since we apply it to the unorganized input point cloud. 
5 We can improve the N(Ck) partition if we know a priori 
that the shape identification procedure described above 
has been applied. Using the information generated during 
shape identification, we construct ,a low order Taylor 
polynomial expansion that equals as a point set a shape 

10 index equivalence class. We begin with an estimate of the 
error that we expect if we expand the signed distance 
function (sdf) in a 1 st order gradient Taylor expansion 
plus quadratic remainder term. Let s (x, y) represent the 
sdf over an open neighborhood N(xo, yo) in the (x, y) 

15 plane. 

Then the Taylor expansion for s in N(x,, y ( ) is 



25 



s(x,y) = s(x 0 ,y 0 ) + 



J_ 
2! 



I^-Uo . y 0 ) • (* - *o ) + f^<*»>. y 0 ) • (y - y 0 )j J 

^(ax t ay)-(x-x 0 ) 2 +-£^(ax,ay)-[(x-x 0 )(y-y 0 )] + 



3 2 S , . / \2 



; a = a(x, y) and 0 < a < 1. 



dy 

We translate s in the z-direction so that s(x 0/ yo) = 
2 0 0, eliminating the constant term. Now rotate the 

orthonormal basis defined by the tangent plane of the sdf 
at (x 0 , yo) plus the z-axis so that i't coincides with the 
orthonormal basis formed the canonical co-ordinate 
system. 



Notes : 
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a. Rotation in the (x,y) coordinate plane will align 

any pair of orthogonal vectors the coordinate 

« 

plane axes. 

b. Rotating the surface so that sfx,, y 0 ) is a maximum 
causes V s(x 0 , y 0 ) = 0. 

Combining the translation and the rotation, the 
Taylor expansion is 



s(x, y) = 



2! 



d z s . ■ . 2 d 2 s , s 2 

(ax,ccy)-x + " x ~ 



dx 2 



dy 



(ax,ay)-y 4 



; 0 < a < 1. 



The coefficients inside the square brackets are the 
principal curvatures at (ax, ay) . They can be computed 
from the mean and Gaussian curvature at (Ox, ay) . 



15 



20 



25 



H (ax, a y) = ^- 



dx' 



(a x;a y) + 



d 2 



d 2 . 



dy 



(a x 9 a y) 



K {a x,a y) = |— s T {ax 9 a y) ^^^(ax 9 a y) 
ox dy 

Using shape analysis, the true value of every term 
in the Taylor expansion is known. The Taylor disk records 
the first two terms in the expansion, the expansion 
point, the disk's radius of convergence, and the maximum 
error incurred on the convergence disk when a sample 
coordinate is approximated by just the constant and 
linear gradient terms. When an intersection curve is 
established, then the definition of : the- curve is # " 
attached. There is a separate attachment for every 
intersection curve. 



44 



WO 2004/057375 




PCT/GB2003/005395 



The a parameter is interesting. Given the center of 
the Taylor disk, the disk radius, and polar angle / 
increments then the a parameter describes the 
intersection as well as the principal curvatures along 
the intersection curve. (A more descriptive name for this 
parameter is "wiref rame" . ) This parameter determines the 
adequacy of the interpolating function in the 3D grid. In 
other words, the adequacy depends on the curvature of the 
parent surfaces along the intersection path. 

We prepare the topological analysis of a shape-based 
surface by creating a manifold whose charts are point 
sets that correspond to shape index intervals. We call 
this a shape index manifold. We explained above how to 
construct the manifold's charts and how to estimate the 
differential properties of a chart. These tasks make 
sense with no external context, i.e., a background 
framework. We refer to properties that make sense in this • 
self-contained context as "intrinsic". Any property of a 
manifold that is not intrinsic is "extrinsic". An example 
of an extrinsic attribute is the relationship between the 
boundaries of a surface and the boundaries of a shape 
index manifold. 

An intrinsic boundary is one that remains a boundary 
under rigid motion. It is a Cartesian tensor by virtue of 
this definition. This boundary separates regions that are 
well approximated by patches taken from a surface of 
revolution. An extrinsic boundary is a boundary that owes 
its existence to the configuration of background 
surfaces. Should the background surfaces move then the 
0. shape of the intersection and indeed the very existence 
of the intersection can change. The choice of solver 
depends on the dataset assumptions . The Shapes ® library 
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does not provide recognition or interaction services 
pertaining to intrinsic boundaries. We agree that a 
characteristic of a sound approach to point set 
classification is a robust algorithm to- compute the 
5 intersection of two extrinsic boundaries . But. we go one 
step further and say that intersection in a region that 
is devoid of extrinsic structure can be still be quite 
complicated, if the intersection involves a non-empty 
subset of a non-dif f erentiable interface that separates 

10 two shapes. 

We seek. a data structure that conveys a framework 
overview. We want it to dontain the complete boundary 
representation for the . framework plus a synopsis of the 
shape of every framework surface and curve. We envision 

15 using this data structure to reply to browser level 

framework data base queries when the caLller does not want 
to open the framework with the standard geometry services 
toolkit. 

We define a structural synopsis to be a Shapes/GQI 
20 topology graph (also known as a boundary-representation 
or u b-rep") plus a shape index manifbld description of 
every 2D node in the b-rep. 

Preferably, a Reeb graph is used to describe a 
configuration's Morse critical points and homotopic 
25 skeleton.. See, M. Hilaga, Y. Shinagawa, T. Komura, T. 
Kunii, Topology matching for full automatic similarity 
estimation of 3D , SIGGRAPH 2001, pg. 203-212, and Silvia 
Biasotti, Topological techniques for shape understanding , 
http://www.cg.tuwien.ac.at/studentwork/CESCG- 
30 2001/SBiasotti/ . We recall the definition of a Reeb 

graph. Let M be a path-co^inected manifold and let f be a 
real-valued on M. Then the Reeb graph associated with (M, 
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f ) is a set of (M, f ) , R) equivalence classes that is 
defined by the relation (x,f(x)) ~ (y, f (y) ) if and only 
if f(x) = f(y) and x, y are members of f-1 ( (f (x) ) . In 
practice, the nodes and arcs in a Reeb graph are - 
5 determined from continuous sampling of homotopic 

identification of height field contours. Since the height 
field is a Morse function, we obtain information 
regarding each' non- degenerate critical point and the cell 
to which the critical point is attached. The discussion 

10 of 2D parameterization herein summarizes the relevant 
facts from Morse Theory. 

Figure 19a shows the Reeb graph of a standard torus. 
Figure 19b schematically illustrates a 2D cell suspension 
that is induced from the axes and planes of symmetry and 

15 critical point theory. The Reeb graph's nodes in this 

case are critical points on the torus, since the critical 
points are isolated and non-degenerate. The graph 
contains 4 nodes. The bottom and top nodes correspond to 
the height field minimum and maximum. These two points 

20 lie on an isoparametric curve of constant positive i 

Gaussian curvature. The other two nodes correspond to the 
lower and upper saddle points. They too lie on an 
isoparametric curve, but this curve has constant negative 
Gaussian curvature. It is easy to see that the symmetry 

25 group of the torus leaves invariant the orbit of a saddle 
point as well as the orbit of the min/max. We remark that 
observing this invariant behavior is an easy way to 
discover symmetry transformations. 

A Reeb graph describes homotopic equivalence. Reeb's. 

30 representation has no concept of shape. The Reeb graph 
says nothing about Gaussian curvature or mean curvature 
or shape index. There is no mention of the two opposing 
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curves of zero Gaussian curvature. Nor is there a 
description of the shape of the saddle point curve that 
bounds the interior void. Homotopic equivalence does not 
leave invariant 2D regions of constant shape index, so it 
5 is not possible to reason about symmetry orbits under 
Reeb equivalence. 

Neither the Hilaga nor Biasotti references augment 
the Reeb graph substrate with curvature information, 
citing the natural instability of curvature estimation in 

10 noise. We agree that curvature estimation in a noisy 

environment is difficult, but we have above that it makes 
sense to treat a noisy surface as a smooth substrate that 
is contaminated with noise. We think of the 
representation of a noisy surface as a minimization 

15 problem in Lagrangian mechanics. The smooth approximant 
represents the kinetic energy in the decomposition. We 
define the curvature of the original surface to be that 
of the smoothed component of the original surface. The 
Riemann integral of the discrepancy between the smooth 

20 approximant and the original surface is a measure of the 
potential energy in the original surface. Noise 
introduces uncertainty into the curvature estimate. We 
compensate for this uncertainty by working with shape 
index intervals, rather than a point-specific value. We 

25 define a shape synopsis diagram (SSD) to a Shapes /GQ I 

boundary representation where a 2D node is a shape index 
manifold. 

We comment further that Hilaga reports the 
development of a similarity metric for a pair of 
30 triangulated surfaces SI and S2 . Their idea is to create 
a level of detail hierarchy for each surface. The authors 
summarize each level of detail by constructing a Reeb 
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graph of the coarsened surface. Hilaga chooses the Reeb 
function to be the integral of geodesic distance measured 
a,t every vertex of SI and S2 . Hilaga approximates the 
integral using Dijkstra's Algorithm. Hilaga argues that 
5 this Reeb function is superior to a height field because 
the integral is insensitive to orientation. 

We want to compare Hilaga' s method to the method of 
the present invention. This is not easy. Hilaga works 
with triangulated surfaces rather than implicit surfaces.. 
10 Consequently we must .be careful regarding the meaning of 
the term "geodesic path" . On a smooth surface we can 
define a geodesic curve to be the straightest possible 
path or the shortest possible path, since the two 
characterizations coincide. On a triangulated surface 
15 they do not. As is common practice, we select 

* straightest possible" geodesies because an existence 
proof is available for * straightest possible" whereas 
none exists if instead we opt for "shortest possible" 
geodesies . 

20 We restrict attention to a single surfade of 

revolution S, because there exists a theorem that says 
that a parallel on S is a geodesic exactly when the 
tangent vector at any point on the surface is parallel to 
S's axis .of symmetry . The symmetry group of S honors this 

25 constraint, so the geodesic is an orbit under the action 
of the symmetry group. 

Figure 20 is cross section of the torus shown in 
Figure 19a. Figure 20 shows the torus opened along a 
planar surface that joins the saddle point circular orbit 

30 and the mi n /max circular orbit. 

In Figure 20., we have overlaid the symmetry group 
orbits associated with critical points. We measure 
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criticality against the standard height field Morse 
function, because geological data is naturally observed 
in depth. We free ourselves of orientation and pose 
limitations by working with the symmetry group orbits of 
5 critical points rather than the isolated critical points 
( themselves. 

Figure 21 is a schematic of the shape synopsis 
diagram (SSD) of the torus shown in Figure 19a. The SSD 
reports locations in normal position coordinates and uses 
10 a homogeneous transformation to correctly position this 
data in model space. We also note that an SSD takes 
little storage beyond that already needed to instantiate 
the b-rep. We have placed details of the shape index 
manifold and shape index chart provided below. 
15 Since an SSD augments the standard b-rep we attach a 

spatial frame to the b-rep node to reference the SSD. 
Since the top-level node in the SSD contains a database 
identifier, standard navigation methods can be used to 
locate the b-rep given the SSD. We call attention to the 
20 arrow notation in the SSD. We use an arrow to define the 
correspondence of ID chart boundaries to 2D charts. The 
arrow that is attached to the concave chart signifies 
that the chart does not bound a sub-volume along its 
outside. Similar remarks pertain to the convex chart. 
25 The Reeb graph of a cyclide is identical to the Reeb 

graph of a torus, as expected since the standard Reeb 
graph does not consider shape. The difference between a 
cyclide SSD and a torus SSD is the absence of circular 
orbits connecting critical points. That is, all. of " the 
30 height field critical points for a cyclide are isolated 
and are fixed by the cyclide' s back to front reflection 
operator . 
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Figure 22 is a cyclide that is shaded according to 
Gaussian and mean curvature. 

Next we, construct the SSD for a 2-torus. Figure 23a 
is a bi-torus with its associated REEB diagram. Figure 
5 23b is the visual representation of the bi-torus shape 

synopsis diagram. A shape index manifold is a patchwork, 
with each patch taken from a surface of revolution. So it 
makes sense to display the surface^ of revolution, where 
each surface is decorated with the bounding curves that 
10 define the patch selection. On the bi-torus 's SSD there 

are two thin circular arcs plus the outer circle identify 
the two toroidal components in the bi-torus. We attached 
the "slab" label to the connective material between the 
two tori, because we' perceive the top' of the bi-torus to 
15 be flat, i.e., its Gaussian curvature is zero. The slab's 
1 internal boundary joins the upper and lower parabolic 

curves on both tori. We have suppressed the ID orbits and 
critical point assignments to keep the diagram readable. 
The bi-torus is. a genus 2 sphere, so there is 1 minimum 
20 and 1 maximum. Since the tori are tilted, the height 

field does not see the any critical points along either * 
saddle point orbit. Finally, the interior surfaces of the 
slab are concave (hyperbolic) in order to conform to the 
torus 's exterior convex surface. 
25 The differences between an SSD and a Reeb graph are 

very clear. The use of homotopic identification means 
that the Reeb graph cannot distinguish homotopic figures 
that have significantly different, curvature, e.g., a 
convex figure from a convex figure. Therefore shape index 
30 analysis is meaningless in the context of a Reeb graph. 

We conclude that although we can map a Reeb graph into an 
SSD, a Reeb graph' cannot support a shape 
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index manifold. 

According to a preferred embodiment of the 
invention,- an efficient hierarchical surface 
5 representation will now be described. We represent an 
implicit surface's signed distance function in a harrow 
band octree encoding of a regularly spaced Cartesian 
grid. We discuss herein how curvature and shape analysis 
of the implicit, surface simplify the octree as a data 
10 structure. Smoothing. or editing in general is likely to 
change a prior shape analysis , so we will also describe 
herein an adjunct shape index representation that enables 
fast updating of the octree's information archive. 

15 1 . We begin by computing the shape index on the 

implicit surface at the finest resolution in the 
octree. 

2 . We partition the shape index data into equivalence 
classes. 

20 a. We find the finest resolution octants that 

contain the surface's height field critical 
points . 

b. We define the following shape index relation on 
the octant corners that were located in the 

25 previous step. 

c . We say that two octant corners are related when 
that their respective shape index values are 
contained in the same Cantzler sub division of 

• - the fundamental space [.-1.0/ +1.0] and that . 
30 both values are bounded away farom both ends of 

the Cantzler subdivision. 

d. We make a second pass over every equivalence 
class of octree nodes that contains less than a 
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threshold number of samples. Let v be such a 
node and suppose that v is on the frontier of 
shape index equivalence classes %1 and %2 . 

e. If %1 and %2 are approximated by a single 

5 elementary shape instance, then we will assign 

the node v to either %1 or %2 . 

f . If %1 and %2 are not approximated by a single 
elementary shape instance, then we assign v to 
%j if the error in the elementary shape 

10 approximation is less than a threshold 

3 . We introduce the shape index information to the 
coarser levels of the octree. 

a. We say that a low-resolution octant is stable 
15 , if there exists an elementary shape 

representation for all highest resolution level 
nodes that are controlled by the coarse level 
octant. 

b. We record stability status. If a low-resolution 
20 node is stable, then it will use the elementary 

shape's algebraic evaluators for signed 
distance, mean and Gaussian curvature, etc. 

c. We recognize that a stable low-resolution 
octree node vhas the same shape resolution as 

25 Vs finest resolution level components. If we 

want a truly coarser approximation, then we 
develop a separate family of elementary shapes 
for each resolution level "in the octree. 

30 4. We prune the stable regions of the octree. 
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a. Given an octant Q at level L such that its leaf 
level resolution is part of an ideal elementary 
shape. 

b. We set a status bit in Q that says it coarsely 
5 resolves part of an ideal shape. 

c. Prune the branch of the octree rooted at Q 

d. If Q does not coarsely resolve part of an 
ideal, then we keep the branch rooted at Q. 
Furthermore, we may need to use the leaf -level 

10* interpolator to estimate curvature. 

The situation most favorable to' this algorithm is 
when the root '' surf ace S* contains a large stable region. 
This can happen if the leaf-level surface S contains a 
15 large region of zero Gaussian curvature, in which case 
the shape index for the region is % or -y 2 . Another 
favorable situation is that both the mean curvature . and 
Gaussian curvature are constant, e.g., a sphere, so that 
the shape index is again a constant. . 
20 Figure 24 is a diagram of the octree with a coarse 

level and leaf level shape index relationship indicated. 
Consistent texture- coding between corresponding regions 
indicates a stable shape index region. 

We consider two examples. The Figure 25 is a diagram 
25 shows part of the French model, A single octree leaf 

defines each plateau. The hemispherical depression has 
constant mean curvature and constant Gaussian curvature, 
so in both regions only one octree leaf is needed. 

Figure 2 6 is an image of the topography of Crater 
30 Lake, Oregon. Figure 26 shows two large plateaus in the 
Southwest part of the lake basin. There are large flat 
regions of the lakebed, so this algorithm will represent 
these regions economically. The small elevation bumps on 
the lakebed will be enclosed in extruded cylinders. 
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We anticipate that the algorithm will perform well 
on the large trimmed. conical plateau in the rear. 

According to a preferred embodiment of the 
5 invention, an implicit surface shape identification 

technique will be described. In particular we show how 
to uniformly approximate an implicit surface by a 
patchwork of smooth shapes, where each shape is a section 
■ . of a .surface of revolution. The implicit surface is not 
10 required to be smooth. We say that the volume between the 
given implicit surface and the patchwork of smooth shapes 
measures the misfit of the approximation. We improve the 
uniform approximation by reducing the misfit. We reduce 
the misfit by approximating the volume as a Riemann sum 
15 * of generalized prisms. 

Shape identification provides a much higher density 
of information. Figure 27 shows a salt weld in the Gulf 
of Mexico. See, M. Hodgkins, M. O'Brien, Salt sill 
deformation and its implications for subsalt exploration , 
20 The Leading Edge, August 1994, pg. 849-851. We enumerate 
the shapes that collectively represent this complicated 
geobody . t 
The shape legend for this image is as follows. 
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Paraboloid cap 
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Toroidal sector 
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Cylinder 
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Barrel 


3 


Toroidal sector 
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Conical sector 
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Barrel 
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Toroidal sector 
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Conical sector 


10 


Cone 



Now we specify the shape identification algorithm. 
Notation 

30 S A noisy implicit surface. 
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S* A smoothed version of S. 

G A 3D regularly spaced grid 

C A grid cell in G with corners {c0 ... c 7 } . 

G The signed distance function defined on C such 

a~ & (0) = S Pi C * 0. . 
•x Tri-linear interpolation associated with grid 

cell C that approximates the signed distance 

function a on C. 
a* The symmetry group of a. 

Shape identification algorithm 

1. / We are given a 3D rectangular or rectilinear grid G 

that contains an implicit surface S. 

2. Construct a smooth surface S* from S, eliminating 
noise, but not low frequency intrinsic shape. 

3 . Sample S* such that there exists a sample at each 
entry and exit point on each grid cell face in G. 

• We know how to do this, since we know the nature of 
a smooth shape that tri-linear interpolation can 
represent . 

4. For each S* sample point estimate the mean and 
Gaussian curvature . 

5. Compute the shape index at every S* sample point. 

6. Mark those samples that are close to the center of 
a Cantzler shape index interval. 

7 . Mark those grid cell samples that are height field 
critical points, i.e., minima, maxima, and saddle 
points. (We know that these occur at grid cell 
corners only.) 

We identify shapes on S*. 

8. For each S* critical point select grid cells 
attached to the critical point that have consistent 
mean curvature, Gaussian curvature, and shape index 
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estimators at all sample points in the grid cell 
These grid cell sets define the initial equivalence 
classes . 

9. We want to define a shape from the region of S* 

5 inside the grid cell equivalence class . We check 

that each admitted grid cell's tri-linear 
interpolant has the same symmetry group. If there 
is no common symmetry group, then we subdivide the 
equivalence class into subclasses that do have a 
10 shared symmetry group. 

10. Using the shape index data and the symmetry group, 
identify a surface of revolution that best fits the 
curvature and unit normal vector field data. 

15 Let C be a grid cell that is unassigned to a grid 

cell equivalence class, but is adjacent to a grid cell 
that has been assigned to a grid cell equivalence class. 
We want a condition to test for admission to the 
equivalence class . 

20 

11. There are two categories of grid cells in an 
equivalence class. A member grid' cell is either 
"shapely" or is a "misfit". 

12. In order for the candidate grid cell to be admitted 
25 as "shapely", it must satisfy the following 

conditions. 

a. Its tri-linear interpolation function has the 
same symmetry group as all other shapely grid 
cell members. 

30 b. If a shapely grid cell in the equivalence class 

shares a grid cell face with the candidate grid 
cell, then they share a smooth boundary across 
the grid cell face, 
c. We re-compute the shape parameter vector and 

35 re-test all grid cells in the equivalence 

class. No "misfit" grid cell is allowed to 
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exceed the misfit tolerance and no shapely grid 
cell can become a misfit. 

13. Test if the candidate grid cell is shapely. 

14. If not, then test if the candidate can be 
considered a "misfit 17 . 

a. A misfit must be less than some threshold and 
no point on the surface sample in the candidate 
grid cell can be further from the ideal shape 
than some other threshold. 

15. If the misfit is acceptably small, then we mark the 
grid cell as a "misfit" and admit the grid cell. 

We specify a misfit reduction algorithm later in 
this part of the description; 

It may be that some grid cells have not been 
assigned to an equivalence class. These grid cells are 
not associated with a shape that has a critical point in 
the smoothed implicit. surface S*. 

16. Form a new equivalence class by. finding connected 
sets of grid cells where the individual grid cell 
tri-linear interpolation functions have the same 
symmetry group and whose shape index samples lie in 
the same Cantzler shape index interval. 

It may be that some grid cells can be assigned to 
more than one equivalence class. If we have a choice, 
then we assign the candidate to the class that accepts 
the candidate as shapely. 

It may be that some grid cells have not been 
unassigned- to an existing equivalence class. Since S* is 
smooth, we know that noise is not a contributing factor. 
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